% 10/01/12 Vasyl Mykahlchuk
%

function [ ] = buildScalarFld( anmFileName )
% Modify matlab path
addpath(genpath('util'));
addpath(genpath('io'));
addpath(genpath('strain'));
addpath(genpath('export'));

% Load animated mesh
anmFileNamePath = strcat('../res/',anmFileName);
anmFileNamePath = strcat(anmFileNamePath, '.anm');
[geometryInfo, topologyInfo] = loadANM(anmFileNamePath);

%{
dimensionsGeometry = size(geometryInfo);
numberOfVertices = int32(dimensionsGeometry(3));
vertexList = geometryInfo(1,1,1:numberOfVertices);

prepareMeshForExport(vertexList, topologyInfo,0);
return;
%}


%==========================================================================
dimensionsGeometry = size(geometryInfo);
dimensionsTopology = size(topologyInfo);

numberOfFaces = int32(dimensionsTopology(1));
numberOfFrames = int32(dimensionsGeometry(1));
numberOfVertices = int32(dimensionsGeometry(3) / 3);
%==========================================================================
% compute adjacency lists
adjList = createAdjacencyList(topologyInfo, numberOfVertices, numberOfFaces);
adjVertList = createVertexAdjacencyList(topologyInfo, numberOfVertices, numberOfFaces);
%==========================================================================

% time-varying scalar field (vertex based)
scalarFields = zeros(numberOfFrames, numberOfVertices, 'double');
% time-varying scalar field (triangle based)
scalarFieldsT = zeros(numberOfFrames, numberOfFaces, 'double');

numberOfFrames = 2;

%
% MEASURE TIME for COMPUTATIONS
%
tStart = tic;

% compute triangle based scalar field
for iFrame = 1:numberOfFrames
    for iFace = 1:numberOfFaces
        tr = topologyInfo(iFace,:);
        
        v1 = geometryInfo(iFrame, 1, tr(1)*3 + 1: tr(1)*3 + 3);
        v2 = geometryInfo(iFrame, 1, tr(2)*3 + 1: tr(2)*3 + 3);
        v3 = geometryInfo(iFrame, 1, tr(3)*3 + 1: tr(3)*3 + 3);
    
        vv1 = geometryInfo(iFrame + 1, 1, tr(1)*3 + 1 : tr(1)*3 + 3);
        vv2 = geometryInfo(iFrame + 1, 1, tr(2)*3 + 1 : tr(2)*3 + 3);
        vv3 = geometryInfo(iFrame + 1, 1, tr(3)*3 + 1 : tr(3)*3 + 3);
    
        %flatten triangles 
        [v1f, v2f, v3f] = fitTrianglesToOnePlane(v1, v2, v3);
        [vv1f, vv2f, vv3f] = fitTrianglesToOnePlane(vv1, vv2, vv3);
    
        strains = tri2D_strain(v1f, v2f, v3f, vv1f, vv2f, vv3f);
        scalarFieldsT(iFrame, iFace) = strains(1);
    end
end

% Elapsed time to compute triangle based deformation scalar field
tElapsed = toc(tStart);

%save('output/log.txt', 'tElapsed', '-ASCII');

tStart = tic;
% compute vertex-based time-varying scalar field
for iFrame = 1:numberOfFrames   
    deformationMax = 0.0;
    for iVertex = 1:numberOfVertices
        i = int32(1);
        scalarValue = 0.0;
        adjTriangle = adjList(iVertex,i);
        while (adjTriangle ~= 0)
            tr = topologyInfo(adjTriangle,:);
            v1 = geometryInfo(iFrame, 1, tr(1)*3 + 1 : tr(1)*3 + 3);
            v2 = geometryInfo(iFrame, 1, tr(2)*3 + 1 : tr(2)*3 + 3);
            v3 = geometryInfo(iFrame, 1, tr(3)*3 + 1 : tr(3)*3 + 3);
            trArea = triangleArea(v1,v2,v3);
            scalarValue = scalarValue + trArea * scalarFieldsT(iFrame, adjTriangle);
            i = i + 1;
            adjTriangle = adjList(iVertex,i);
        end
        if scalarValue > deformationMax
            deformationMax = scalarValue;
        end%if
        scalarFields(iFrame, iVertex) = scalarValue;
    end%for
    for i=1:numberOfVertices
        if deformationMax ~= 0
            scalarFields(iFrame, i) = scalarFields(iFrame, i) / deformationMax;
        end
    end%for
end

%
tElapsed = toc(tStart);
%save('output/scalarFields.txt', 'scalarFields', '-ASCII');

% compute average vertex-based scalar field from time-varying scalar field
avrgScalarField = zeros(numberOfVertices, 1, 'double');
for iVertex = 1:numberOfVertices
    avrg = 0.0;
    for iFrame = 1:numberOfFrames
        avrg = avrg + scalarFields(iFrame, iVertex);
    end
    avrg = avrg / numberOfFrames;
    avrgScalarField(iVertex) = avrg; 
end
%save('output/avrgScalarField.txt', 'avrgScalarField', '-ASCII');

targetScalarField = scalarFields(1,:);
[sizeX, sizeY] = size(targetScalarField);
targetScalarField = reshape(targetScalarField, sizeY, sizeX);
%targetScalarField = avrgScalarField;
criticalPoints = extractCriticalPoints(targetScalarField, adjVertList);
%save('output/criticalPoints.txt', 'criticalPoints', '-ASCII');
colorList = createColorInfo(criticalPoints, targetScalarField);
vertexList = geometryInfo(1,1,1:numberOfVertices*3);
prepareMeshForExport(vertexList, topologyInfo,colorList,anmFileName);

end%function buildScalarFld

%========================================================================
% Create color info based on of deformation scalar field
function [colorInfo] = createColorInfo(criticalPoints, scalarField)
[nV, x] = size(scalarField);
colorInfo = zeros(nV, 3, 'uint8');
for i=1:nV
    colorInfo(i,:)=[255 255 255];
end%for
[nCP,x]=size(criticalPoints);
for i=1:nCP
    color = uint8( 2550 * scalarField(criticalPoints(i)) );
    if color > 255
        color = 255;
    end
    colorInfo(criticalPoints(i),:)= [100 255 color];
end%for

end%function create color info

%========================================================================

function [vertexList,faceList,vertexColorList ] = prepareMeshForExport(geometryInfo, topologyInfo, colorInfo,anmFileName)
dim = size(topologyInfo);
nF = dim(1);
[dim1,dim2,nV] = size(geometryInfo);
nV=nV/3;

%If color info is empty, then construct default one with (255,255,255)
if colorInfo == 0
    colorInfo = zeros(nV,3);
    for i=1:nV
        rA = rand(1:3);
        colorInfo(i,1) = int32(rA(1)*255);
        colorInfo(i,2) = int32(rA(2)*255);
        colorInfo(i,3) = int32(rA(3)*255);
    end%for
end%if
vertexColorList=colorInfo;

%Construct face list
faceList = zeros(nF, 3);
for i = 1:nF
    faceList(i,1)=3;
    for j=1:3
        faceList(i,j+1)=topologyInfo(i,j);
    end
end%for

%Construct vertexList
vertexList=zeros(nV,3);
nV = nV - 1;
for i=0:nV
    for j=1:3
        vertexList(i+1,j) = geometryInfo(1,1,3*i+j);
    end%for
end%for

anmFileName = strcat('../target/', anmFileName);
exportFilePath = strcat(anmFileName, datestr(now, '_mm-dd-yy_HH-MM'));
exportFilePath = strcat(exportFilePath, '.ply');
exportPLY(exportFilePath, vertexList, vertexColorList, faceList);

end%function prepareMeshForExport